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We address the existence and properties of discrete embedded solitons (ESs), that is localized 
waves existing inside the phonon band in a nonlinear dynamical-lattice model. The model describes 
a one-dimensional array of optical waveguides with both x' 2 ' (second-harmonic generation) and 
x' 3 ' (Kerr) nonlinearities, for which a rich family of ESs are known to occur in the continuum 
limit. First, a simple motivating problem is considered, in which the nonlinearity acts in a 
single waveguide. An explicit solution is constructed asymptotically in the large-wavenumber limit. 
The general problem is then shown to be equivalent to the existence of a homoclinic orbit in a 
four-dimensional reversible map. From properties of such maps, it is shown that (unlike ordinary 
gap solitons), discrete ESs have the same codimension as their continuum counterparts. A specific 
numerical method is developed to compute homoclinic solutions of the map, that are symmetric 
under a specific reversing transformation. Existence is then studied in the full parameter space 
of the problem. Numerical results agree with the asymptotic results in the appropriate limit and 
suggest that the discrete ESs may be semi-stable as in the continuous case. 
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I. INTRODUCTION 



An embedded soliton (ES) is an isolated solitary wave in a non-integrable system that resides insides the 
continuous spectrum of linear waves. Unlike regular gap solitons, the existence of ESs in continuum models is 
generally of codimension-one in the model's parameter space. That is, they are embedded into wider families of 
generalized solitary waves which have nonvanishing periodic tails attached to them, such that the ES occurs at 
isolated values of the frequency, or other intrinsic parameter of the solution family at which the tail's amplitude 
vanishes. The concept of ESs was introduced in Q (see also 0,0 for more details), although particular examples 
of waves of this kind were actually known earlier (but not under this name) 0,0,0- In the last few years, 
ESs have been shown to be of relevance in a diverse range of nonlinear- wave models, see, e.g., 7, 8, 9, 10] and 
references therein. One of their characteristic properties is that they arc typically, at best, semi-stable. That is, 
the linearization around them gives rise to no unstable eigenvalues, but a non-trivial zero-eigenvalue mode can 
be found. The latter leads to an algebraic-in-time instability for energy-decreasing perturbations, and a similar 
algebraic-in-time relaxation back to the original pulse for energy-increasing perturbations . Only in special 
cases of systems which have more than one dynamical invariant, or for which the travelling- wave svstem is 
reversible but not Hamiltonian, have examples been found of ESs that are truly asymptotically stable [l0|. Also 
known are physically relevant models with special symmetries that support continuous (rather than discrete) 
families of ESs, and a large part thereof may be stable [TlLIT^ . 

All the above-cited examples pertain to continuum models. A fundamental issue is whether solitons of the 
embedded type may also occur in discrete nonlinear equations modelling dynamical lattices. For example, in 
PH it was argued that moving topological defects in Frenkel-Kontorova (discrete sine-Gordon) lattices can 
meaningfully be described as embedded kinks. Some of specific examples of moving discrete solitons constructed 
by means of the inverse method (a model is sought for to which a given moving soliton is an exact solution [l5^ 
may also have the embedded character. Such waves, which connect a zero state to itself, translated through 
a multiple of the period of the lattice potential (the 'topological charge'), were first (indirectly) identified by 
Peyrard and Kruskal ^3] an d have later received much attention in a wide variety of applications [is|. In 
a co-moving frame, these waves should be thought of as localized solutions to infinite-dimensional advance 
and delay (nonlocal) equations, see works 0, |2(j where branches of embedded kinks with the topological 
charge 1 and 2 were computed. Another recent study [2l| reported the existence of ESs in a discrete lattice 
of the Ablowitz-Ladik type, but with a quintic nonlinearity added to the usual cubic terms, and also with an 
additional next-nearest-neighbor linear coupling. By tuning the nonlinear terms, an explicit analytical solution 
for a discrete ES was found in that model (which somehow resembles the above-mentioned inverse method). 
Apart from these, we are aware of no other work that systematically considers the issue of ESs in discrete 
models. 

A general objective of this paper is to understand the existence and multiplicity properties of discrete solitons 
in lattice systems for which the corresponding continuum model is known to support ESs. The first fundamental 
issue is whether the codimension of their existence is preserved under discretization. In dissipative systems, 
solitary pulses are supported due to the balance between forcing and damping, hence the corresponding ho- 
moclinic orbits are always of codimension one. Thus, formally similar to ESs, dissipative solitons typically 
exist as solutions to the corresponding stationary ordinary differential equation at discrete values of a relevant 
parameter (temporal or spatial frequency, depending on the physical realization). Then, it will generically hap- 
pen that, under the discretization of that equation, the codimension of the solitons decreases. This is possible 
because discretization of homoclinic orbits to hyperbolic equilibria leads to transverse homoclinic connections 
that exist for a range of parameter values j2^|. However, for ESs we shall show that, in contrast, discrete ESs 
typically keep the same codimension as their continuum counterparts, i.e., they also exist at discrete values of 
the corresponding physical parameter. 

In this work, we restrict ourselves to a discretization of a single continuum model, for which the discrete 
counterpart finds a direct physical interpretation. In fact, our results and methodology will carry over to other 
discretized models that support ESs (this expectation is supported by the general fact that equations of the 
discrete nonlinear-Schrodinger type , that we consider below, may be derived as an asymptotic limit from a much 
larger class of discrete models [l6(). Specifically, we shall start with the continuum model in which ESs were 
identified [p. It describes an optical medium with competing quadratic (x' 2 ') and cubic (x'" 3 '') nonlinearities 
(see [23l l24l |25| and references therein for physical applications). In this model, the evolution variable z is the 
propagation distance in the nonlinear optical waveguide, and the variable t is either the reduced time (in the case 
of temporal solitons) or, which is more physically relevant, a transverse spatial coordinate in a planar waveguide. 
In the temporal case, such a model can be written in dimensionless form as (see detailed explanations in reviews 
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where M and v are the local amplitudes of the fundamental-frequency (FF) and second-harmonic (SH) fields, an 
asterisk denotes complex conjugation, 6 is the relative dispersion of the SH, q is the wavenumber mismatch, and 
7i,2 measure the relative strength of the cubic (Kerr) nonlincarity compared to its quadratic (SH-generating) 
counterpart, whose strength is normalized to be 1. In the physical model, 71.2 must have the same sign 
(typically, they are positive, but may, in principle, be negative too). The spatial variant of the system takes 
the same form, with a p artly different interpretation of the dimensionless parameters, and is easier to implement 
experimentally |23l l24l |25| . In the spatial domain, one always has 5=1. 

The existence of ESs in this model was established in Q,0| for both cases, 712 > and 71.2 < 0. In the case 
of relevance to spatial waveguides (5 = 1), fundamental ESs are absent in the model, while higher-order ones 
can be found. 

In models such as |QJ which includes x^ interaction, ESs can be embedded only in the continuous spectrum 
of the SH component, while the FF wavenumber can never be located inside the continuous spectrum. This 
feature is stipulated by the asymmetry between the two equations in system (Q. If the SH supports linear 
waves, while the FF has the possibility of exponential localization like e _A '*', then the u 2 term, which drives 
the SH field in the second equation JTJ, allows it to have tails that decay exponentially at a rate e" 2A '*L Such a 
solitary wave was said in reference [l| to be tail-locked and, accordingly, the SH equation to be nonlinearisable 
(because the quadratic term can never be neglected in this equation). The tail locking does not take place if 
the SH field supports its own exponentially decaying tails that asymptotically (for \t\ — > 00) dominate over the 
tails induced by the term u 2 . 

In this paper we study the existence of discrete ESs in a model arising from discretization of the (actually, 
spatial) i-coordinate in One motivation for doing this is to understand the effect of numerical approximation 
(which also implies discretization, once a finite-difference scheme is employed) on the computation of ESs in 
model systems. But this is not our primary motivation. Lattice models play an increasingly important role 
in describing physical phenomena in a number of newly developed experimental settings. In particular, the 
discrete version of the x^ '■ model describes an array of the corresponding linearly-coupled waveguides. 
The creation of discrete spatial solitons in a system of channel waveguides with the x nonlinearity was recently 
reported |27|, Our assumption concerning relevant solutions is that, once the continuum version of the model 
readily gives rise to ESs, then there is a chance to find them in the corresponding lattice model too. In this 
paper, we show that this is the case indeed and investigate how we may then pass to the continuum limit. We 
report the corresponding discrete ES that can be found, under special conditions, in an approximate analytical 
form, and, in the general case, numerically. 

The paper is organized as follows. In Section 2, we present the lattice model to be considered in this work and 
discuss its physical applicability. In Section 3, an asymptotic analysis is undertaken for a reduced model where 
the cubic nonlinearity is present only at a single site. Section 4 then goes on to discuss the meaning of discrete 
ESs, as a solution to a system of stationary finite-difference equations, in terms of a homoclinic orbit to a fixed 
point of a discrete-time reversible map. This makes it possible to understand that the solutions we seek must 
have codimension one, similar to ESs in continuum models. In terms of the same approach, in Section 5 we 
describe the numerical procedure developed to search for discrete ES solutions. Note that this requires special 
adaptation or other methods for finding homoclinic orbits of maps, owing to the non-hyperbolic nature of the 
fixed point. Our approach heavily relies on the reversibility although it can be readily modified to treat non- 
reversible ESs. Numerical results are reported in Section 6, in the form of the corresponding codimension-one 
solution branches in the parameter space. An array of actual shapes of the ESs is displayed too. Our numerical 
results also show that the asymptotic analysis of Section 3 can well predict the codimension-one set of ESs in 
the parameter space in the model when the wavenumber is large. The paper is concluded by Section 7 along 
with a preliminary stability result for the fundamental ES. 
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II. THE MODELS 



We now consider Eqs. Q in the spatial domain. Direct discretization of the spatial second derivative produces 
a lattice model 



+ ]-D(u n+ i + u n -i - 2u n ) + u*v n + 7i(|u„| 2 + 2\v n \ 2 )u n = 0, 
az A 

i ~^z + \ 5D ^ Vn+1 + Vn ' 1 ~ 2w ") + qVn + + 2 T2(b«| 2 + 2\u n \ 2 )v n = 0, 



(2) 



where n is the discrete transverse coordinate which assumes integer values, and D = 1/h 2 with h the stepsize 
of the discretization. Effective coefficients of the lattice diffraction for the FF and SH waves are Z?i = D and 
£>2 = SD; with 8 = 1, this yields D2 = D\ in Eqs. J2J. Unequal values of these coefficients, D\ ^ D2 (i.e., 
5 =/= 1), including those with opposite signs, can in principle be realized experimentally (the latter possibility was 
recently employed to predict the existence of discrete gap solitons 28] ). To do this one would use a diffraction- 
management technique, which is based on oblique propagation of the beam across the waveguide array which 
is modeled by the discrete system (2^, [3(j • In that case, the effective lattice-diffraction coefficients are 

D m = D<£> cos(mQ), m = 1 , 2, (3) 

where Q = fcsinf is the transverse wavenumber, k is the beam's wavenumber (normalization is such that the 
array spacing is unity), is the angle between the Poynting vector and the coordinate z running along the 
waveguides, and ufX are the original lattice-dispersion coefficients, corresponding to Q = 0. As seen in Eq. 
one can efficiently control the size and signs of D% t 2 by choosing an appropriate angle 9. In particular, the ratio 
D 2 /Di can be made large by choosing Q = tt/2 + e, or Q = 3n/2 — e, for a small positive e. 

To cast the model in a normalized form, we note that D± can be made positive, if it was originally negative, 
by means of the complex conjugation (i.e., Eqs. @ are replaced by their counterparts for u* L and v n ), combined 
with the changes v n —* —v n and 71 2 — * —71.2- The size of D\ may be set equal to 1 by means of the rescaling: 
zD\ —* z, (u n , v n ) / D\ — > («„, v n ), which leads to the final normalized form of the discrete model, 



i ~([z + \^ Un+1 + Un ~ 1 ~ 2 '"™' ) + U * nVn + ^d"""! 2 + 2 \ v n\ 2 )u n = 0, 

Cj/tJnrt 1 _ , v 1q / I I Q I I ^ \ 

h »o{v n +i + Vn-1 - 2v n ) + qv n + -U n + 2j 2 {\Vn\ + 2 \U n \ )v n = 0, 

az 2 2 



(4) 



where 



q = q/Di, 71,2 = £>i7i,2- (5) 

The dimensionless constants (jSJ, including 8, may, in general, be positive, zero, or negative. 

To obtain some analytical results, we shall also introduce a simplified model where ESs can be found in 
approximate form. In the simplified model, only the central site (the one corresponding to n — 0) carries the 
Kerr nonlinearity. Such a model is physically feasible too, if the central waveguide in the array is doped to 
enhance its Kerr nonlinearity. Thus, the simplified model is based on the equations 

i ~J 11 + \ ( u ' n+1 + Un ~ 1 ~ 2u ") + u * nVn + ^ l£ " (' M °' 2 + 2 I W °| 2 ) u ° = °> 

(6) 

i ~dz + 2 " +1 + Vn ~ 1 ~ 2 + ^ Vn + 2 U " + 2 ^ 2£ " ' + 2 ' V ° = °' 

where e„ = for n ^ 0, and e n = 1 for n = 0. 

In either model, or (JSJ, the ES corresponds to a solution with FF component exponentially localized: 

u n ~ ^4exp(ifcz — A|n|) as \n\ — * oo, (7) 

where A is a real amplitude and A is positive and real. Simultaneously, the propagation constant k must belong 
to the phonon band of the SH equation. That is, the SH component of the solution generically will have a 
nonvanishing tail, of the form 

v n ~ C exp(2ifcz — ip\n\) as n| — > oo, (8) 
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where p is real. The existence of an ES corresponds to the vanishing of the constant C in Eq. © . Linearization 
of Eqs. © an d @ an d subsequent substitution of expressions © and © yield the following relations: 



k = 2sinh 2 (A/2), (9) 

g- 2fc = 2<5sin 2 (p/2). (10) 

From Eqs. © and (|10|) we see that a discrete ES may exist in the case of k > 0, with g belonging to the region 

0<q-2k<2S or 2S<q-2k<0, (11) 

depending on whether <5 is positive or not. Note that we can only choose k, q and d independently in the set of 
(k, q, 5, \,p) since A and p are determined by the other three parameters via Eqs. © and <|1U|) . 

III. ASYMPTOTIC ANALYSIS FOR THE SIMPLIFIED MODEL 

We first consider the simplified model ©, with the nonlinearity present only at the central site. An ES 
solution is sought asymptotically under the assumption that \v n \ -C \u n \ = 0(1), the validity of which will be 
checked a posteriori. Accordingly, the quadratic term in the first equation of the system © may be dropped to 
first approximation, which makes the expressions © and © an exact solution for n^O. At the point n = 0, 
the cross-phase-modulation term, |^o| 2 uo; in the first equation © may be dropped too to leading order. Then, 
the equation at n = yields a final result for the FF component of the soliton, A 2 = 7-f 1 sinh A, which implies 
that the discrete soliton in the FF component is supported by itself, without coupling to the SH component. 
Further, using Eq. © one can express A 2 in terms of k, 



i 2 = 7 rV^ + 2), (12) 

which means that the solution exists only in the case of 71 > 0. 

Now, we tackle the second equation of © and substitute the FF field in the form of Eqs. © and ©. 
Obviously, at n ^ 0, an exact solution can be found, which precisely corresponds to an ES, as it does not contain 
the nonvanishing tail © and is 'tail-locked' to the square of the FF field, as explained in the Introduction. We 
find explicitly that 

v n = Bexp{2ikz - 2X\n\), (13) 

where 

B. ^ ^ s -^ J^L^. ,14) 

2 [2(5 sinh 2 X+(q- 2k)} 2 7i 28 ■ k(k + 2) + (q- 2k) 

Substituting the expressions © and l|13f) into the second equation of © at n = 0, we have 



coshA^I? 2 -. (15) 
d 71 

Here, following the assumption \v n \ -C \u n \, we neglected the self-phase-modulation term |wo| 2 vo at this point. 
Equation l|15fl also implies that 

P- > 1- (16) 
071 

and, especially, that 5 has the same sign as 72 since A, 71 > 0. Finally, using the relations © and i|15|) . we 
obtain 

72 
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which is positive by virtue of (|16f) . This result selects the single value of k at which the simplified model admits 
the existence of an ES, with tail- locked SH component. It follows from Eq. (|17|l that such a solution exists if 
the relative lattice-diffraction coefficient (see Eq. |[SJ) belongs to the interval 

< \5\ < 2 7l /| 72 |, 71 >0, (18) 

where S and 72 must be of the same sign. 

It is now necessary to check compatibility of the solution with the underlying assumption, \v n \ -C |tf n |, which 
implies A 2 <C B 2 (see Eqs. (|12[1 and l|14|)). Straightforward consideration shows that this condition amounts to 
1^1 *C |'Ti,2 1 - We also note that, if the FF component of the soliton is strongly localized, so that e~ A <C 1 (see 
Eq. 10), the simplified model is actually equivalent to the full one, as the nonlinearity in the first equation of 
J2J is then negligible at n 5^ 0. From JjJJ we see that this condition implies that k must be large and hence, 
from JTJJ, tnat \ S \ < 1 - 

A noteworthy feature of Eqs. I|17|) and (|18|l is that they do not involve the mismatch parameter q. However, 
the condition that the soliton found above is embedded implies that k given by Eq. (|17|l must belong to the 
interval (| 1 1 p . which relates k to q, as \S\ <C 1. This approximate analysis, valid too for the original model in the 
limit of small \6\, means that, for an isolated selected k- value given by l|17|). there exists a curve of single- humped 
ESs approximately parametrized by q in the narrow interval 

q G (2k, 2k + 28) for S > or (2Jfc + 25, 2Jfe) for S < (19) 

in the (q, /c)-plane for <5 fixed. 

IV. REVERSIBLE MAPS 

Let us now consider the model @ in the general case, and understand, from a dynamical systems point of 
view, why discrete ESs should exist. In this context, one of essential issues is the transition to the continuum 
limit. The scaling leading to Eq. in which D\ = 1, does not allow this. Hence, we undo this scaling and 
set Di = D, Z?2 = SD. As above, we look for stationary solutions in the form 

u n - U n e lk \ v n = V n e 2lkz , (20) 

where U n and V n are real. Equations J2J thus reduce to 

1 



^D(U n+1 + LT„_ a - 2U n ) - kU n + U n V n + 71 (U 2 + 2V 2 )U n = 0, 
\5D(V n+1 + F n _! - 2V n ) + (q- 2k)V n + ~U 2 + 2 l2 (V 2 + 2U 2 )V n = 0, 



(21) 



where parameters without the tildes are related to those with tildes by 

q = qD, k = kD, 7l , 2 = % 2 /D (22) 

Note that the possible existence regions of ESs, given by expression Qll|). now becomes 

< q - 2k < 2SD or 25D <q-2k<0. (23) 

First and foremost, we want to establish the codimension of ES solutions to l|21l) . In order to do that, it is useful 
to recast the system as a four-dimensional map. Specifically, upon scaling the variables as £ n = y/2\^i\/D U n 
and r\ n = y/2\ji\/D V n , we can indeed view Eqs. 1)2 l|l as defining a four- dimensional map, 

Cn+l = F e (C„), (24) 

where ( n = (Xn,€n, Mn, Vn) and 



Fe(Cn) 



( Zn \ 
Vn 

V/i 2) (Cn)/ 



G 



Here we define 



/ e (1) (Cn) = -Xn + 2»itn- tenth + Kl{£ + 2T&)£ n ), (25) 
/ e (2) (Cn) = -Mn + 2^„ - J" 1 f ~e£ + 2«a(»£ + 2£> n J , (26) 



where 



k 2k -q 

v\ = 1 + — , ^2 = 1' 



D SD ' V ^iTil 

V (27) 



Ki = sgn(7i), «a=7a/|7i|- 



The map F e is reversible [31132], in the sense that F" 1 (i?(C„)) = i?(F e (C„)) holds, where i? : R 4 -> M 4 is the 
(linear) involution given by 

). (28) 

Note that the map F e is also reversible, under the second reversibility 

R : M n) -> (/ e (1) (Cn),^,/i 2) (Cn),t^). (29) 

In what follows, we consider orbits of the map that are reversible with respect to R. These will correspond to 
ESs that have their central peaks on a lattice site, which are the kind of lattice solitons that are most frequently 
observed in physical applications. In contrast, ESs that are symmetric under R would have their central peak 
between two lattice sites. Both R and R arise from the fact that the ordinary differential equations for steady 
solutions of the continuum model QJ, 



ld 2 U 



2 dt 2 



kU + UV + 71 (U 2 + 2V' 2 )U = 0, 

(30) 



+ iq ~ 2k)v + \ u2 + 2l2{y2 + 2lj2)v = °' 

are reversible too, under an involution R : (U, dll/dt, V, dV/dt) i— > (U, —dU/dt, V, —dV/dt). 

A fundamental characteristic of reversible maps is that if {Cn}S£L-oo 1S an orbit then {i?(C- n )}^L_oo is also 
an orbit. We say that an orbit {Cn}5^L-oo 1S symmetric (with respect to the reversibility) if Cn+i = R(C-n)- 

Denote Fix(i?) = {C £ M 4 F £ (Q = R{()}- We easily see that C = (x, t, A*j v) 6 Fbc(R) if and only if x = /e (1) (C) 
and fi — fe (C)- Thus, Fix(i?) depends on the particular form of the map F e . An orbit {Cri}5J=_oo i s symmetric 
if and only if Co £ Fix(i?), i.e., £_i = £i and 77 1 = 771. 

The map F e has a fixed point at the origin O, whose Jacobian matrix is 



J 



f 1 

-I 2z/i 

1 

V 0-1 2^ 2 



It follows from Eq. (|23|) that v\ > 1 and j^l < 1, hence the origin is a fixed point of F e of saddle-center type. 
The saddle-center has one-dimensional stable and unstable manifolds, W s (0) and W u (0), which are tangent 
to the stable and unstable subspaces spanned by the vectors (1, v\ — \J v\ — 1, 0, 0) and (1, v\ + \/v\ — 1, 0, 0), 
respectively, and a two-dimensional center manifold, W c (0), which is tangent to the center subspace spanned 
by a set of two vectors, (0,0,1,0) and (0,0,0,1). By the fundamental property of reversible maps, W s (0) = 
R(W u (0)) and W u (0) = R(W s (0)). Thus, if there exists an orbit {Cn}£L-oo on W u (0) such that Co £ Fix(i?), 
then it is also contained in W s (0) and is a symmetric homoclinic orbit to O. 

If F e were not reversible, then such intersections between the one-dimensional stable and unstable manifolds 
in the four-dimensional phase space would be of codimension two. However, symmetric homoclinic orbits are 
of codimension one, since, for their existence, we require an intersection between the one-dimensional unstable 
manifold W u (0) and the two-dimensional manifold Fix(i?). Thus, since a homoclinic orbit to O for F e represents 
precisely an ES in the lattice system J5lJ, we have established that: 



7 



Embedded solitons of the lattice model are of codimension-one in the parameter space, provided they are 
symmetric under a reversibility equivalent to R (or R). 

Note that this is identical to the multiplicity result known in the continuum version of the model 0] • 

Finally, in what follows we shall also treat the case of pure x^ nonlinearity, 71 = 72 = 0. This case is 

physically important, as experiments could be quite conceivably be set up in a medium with negligible Kerr 

nonlinearity. However, without the x^ terms, no ES exists in the continuum model 0. It will therefore be 

important to find out whether the same is true in the discrete model too. 

With 71 = 72 = 0, the scaling leading to F c becomes invalid, so we shall consider instead a family of 

four-dimensional maps 



Cn+l — G e ((n), 



(31) 



where we define 



Ge(Cn) 



( U \ 

Vn, 



with 



9P(Cn) = 



- Xn + 2vx£ n - (e£ n ri n + (1 - e)(£ + 2»£)£ n ), 

-/in + 2"2Vn - 5- 1 Qe£ + 2(1 - e)(ril + 2^)^ , 



(32) 
(33) 



with V\ and V2 given by (|27|l . The map G e is reversible under the same involutions R and R. The map G\ is 
equivalent to 121|) with 7 X = 72 = if one sets £ n = U n and rj n = V n , and simultaneously Go coincides with 
Fq with Kx = K2 = 1. The origin O is also a saddle-center of G e and has the same stable, unstable and center 
subspaces as those of F e . So we can apply all the arguments presented above for F c to G c if we replace F e with 
G e in the definition of Fix(i?). In particular, a homoclinic orbit to O for G e represents an ES in the lattice 
system (J2J) with 71 =72 = 0. 



V. NUMERICAL PROCEDURE 

Several numerical procedures exist for finding homoclinic orbits to fixed points in maps, see, e.g., [36ll3^ .l38l 
I39I Eo) . Here we shall use an adaptation of these methods that takes into regard both the reversibility and the 
non- hyperbolic nature of the fixed point. 

To compute symmetric homoclinic orbits for F £ , we consider the three-dimensional algebraic problem 

X-n - (vi - £-jv = 0, (34) 

£l=£-i) Vi=V-i, (35) 

for N > sufficiently large, where C±i = {x±ii C±i> M±i> ^±1) = F^ f±1 {x-N 1 £-n, 0, 0). The condition l|3^|) 
means that the point (x-jv; C-JV; 0j 0) lies in the one-dimensional unstable subspace of the fixed point at the 
origin, and condition l|35f) means that (0 = (xo, £oj Moj Vo) G Fix(iZ), i.e., the finite orbit {Cn}n=-N intersects 
Fix(i?) at n = 0. Thus, adding a parameter as an extra unknown variable, Eqs. (I34|l and i|35|) represent a 
formally well-posed system of three equations for two unknowns, X-n and £-JV, and the free parameter, which 
can be chosen to be any of e, 5, fc, D, q, k\ or K2- A non-trivial solution for N sufficiently large gives an 
approximate homoclinic orbit {Cn}n=-N °f that is symmetric under R, which in turn corresponds to a 
discrete ES solution of the lattice equation @ for 71,2 7^ 0. In the case 71 2 = 0, the same treatment can be 
used to find symmetric homoclinic orbits for G e . Note that one test of the validity of this approximation is to 
measure the distance of the point (x-jv 5 £—N, 0, 0) from the origin. By virtue of the stable manifold theorem 
for maps [3j|, we know that the error will be proportional to this distance squared. 

Branches of solutions to Eqs. i|34|) and (|35|l can be continued in a second parameter using pseudo-arclength 
continuation. To this end, we use the code AUTO j34jj. The problem now is finding a good initial point along a 
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-1.2 1.2 

Cn-1 

FIG. 1: Homoclinic tangles for the two-dimensional map / with v\ = 1.25 and k\ = 1, drawn using the software Dynamics 
|35|. The circle "•" represents the saddle at the origin. 

branch of discrete ESs. Here we can use a regular perturbation approach by first finding solutions to the map 
.Fo, making use of the fact that when e = 0, the {xm £n)"Pl ane is invariant under F e . The restriction of Fq onto 
the invariant plane is given by 

(Xn+l,£n+l) = f(Xn,£n), (36) 

where 

f(Xn,^n)= (£n,-X« + 2l/i£„ - 

The two-dimensional map / also has a hyperbolic saddle at the origin, and is reversible under an involution, 

R ■ (Xn,£n) >-> (£r»,Xn)- 

The stable (resp. unstable) manifold of the saddle, W s (0) (resp. W u (0)), is tangent to its stable (resp. 
unstable) subspace spanned by a vector (l,fi — \/ v\ — V) (resp. (1, t>i + \Jv\ — !))■ By the reversibility, 
W /s (0) = fl(W u (0)) and vice versa. When v\ > and «i > 0, the stable and unstable manifolds intersect 
transversely and form homoclinic tangles, as shown in Fig. ^ 

Let N > be a sufficiently large integer, and let (x-n, £-n) be a point on the unstable subspace such that 
(6v+i,XJV+i) = (%-n,£-n), where {x~n+i,£n+i) = f 2N+1 (x~N, £-n)- By the reversibility of/, the point 
(Xn+i,£,n+i) must be contained in the stable subspace. The two points (x-jsr,^-N) an( i (xn+i, £n+i) are also 
close to the saddle when N > is large. Hence, the orbit leaving at (x-n,£-n) on the unstable subspace and 
arriving at (xn+i, £,n+i) on the stable subspace is a good approximation to a symmetric homoclinic orbit of /. 
Using an adaptation of the driver HomMap [sij l4lj to AUTD97, we can easily find such approximate homoclinic 
points. 

Now, in order to find non-trivial symmetric homoclinic solutions of F e for e > 0, we take e as the additional free 
parameter and choose (x-n, £,-n, e) = (X-JVj £-jv, 0) as the starting solution to H^4|l and where (x-n, C-Jv) 
denotes the homoclinic point on the unstable subspace for /, obtained using the above procedure. Fixing m = 1, 
we then performed continuation of these algebraic equations in AUTO, using 5 as the continuation parameter. To 
obtain symmetric homoclinic orbits for n\ = —1, we also take 5 and K\ as the free and continuation parameters, 
respectively, and continue the solution obtained above for m — 1 and e = y2/ (-D|7i|). The results are presented 
in Figs. 12141 As shown in Figs. |2 and |3 new branches bifurcate from the one with e = at discrete values of 
5 and can be continued to e = y/2/ (D[yi\) by varying S and e. As shown in Fig. 0] the branches were also 
continued from m = 1 to Ki = — 1 by varying S and k±. For k\ = — 1 and K2 = 1, we could not find such 



9 



0.5 - 



w 



-0.5 



. (a) 










I , i , , , , i , , , , i 



w 



1 2 
5 




FIG. 2: Numerical continuation of symmetric homoclinic orbits for F £ when D — 100 and N = 85: (a) k = 0.3, q — 5 
and Ki — K2 = 1; (b) k = 0.5, g = 1, k\ = 1 and «2 = — 1. Here e and 5 are varied. The dotted line represents 
e = ^2/ (DM) with 71 = 0.05. 
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FIG. 3: Numerical continuation of symmetric homoclinic orbits for G e 
and S are varied. The dotted line represents e = 1. 



when D = 10, fc = 3, q = 1 and iV = 15. Here e 



symmetric homoclinic orbits of F e with sufficient precision. In these computations, we also chose the value of 
N such that the distance between the approximate homoclinic point and the origin was 1.5 x 10 -3 at most, and 
checked that the results did not change significantly under increase of N. Our computations also suggested that 
a possibly infinite number of branches of symmetric homoclinic orbits could be obtained when k or q is large or 
S is small. These branches show oscillations in the parameter plane, which are sensitive to the value of N and 
do not appear to converge as N — -> 00. These, probably spurious branches, result from a large ratio between 
the imaginary eigenvalue of the linearization and the real eigenvalue, which is well known to be a singular limit 
for equations that bear ES 0>E3> an d needs to be treated with great caution. In the results that follow, we 
have stopped computation at points where such oscillations first become evident, and have checked all results 
for consistency in the limit of iV — > 00. 



VI. NUMERICAL RESULTS 



We now present continuation results obtained with AUTO for symmetric homoclinic orbits of F t (or G e ) under 
simultaneous variation of two relevant parameters within the saddle-center parameter regions. By varying the 
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FIG. 4: Numerical continuation of symmetric homoclinic orbits for F e when D — 100, k = 0.5, q = 1, «2 = — 1, N = 85 
and e = 0.6324555 ft! \/2/ (DpyiJ) with 71 = 0.05. Here S and «i are varied. The dotted line represents Ki = — 1. 

discreteness coefficient D up to large values, we are also able to compare the results to those of the continuum 
limit, D — » 00. We shall treat the cases 5 > and <$ < separately and also discuss the possibility of discrete 
ESs in the absence of cubic nonlinearity. 

A. The case of 8 > 

Figure [5] depicts branches of discrete ESs in the presence of terms, with 71 = 72 = 0.05 and k — 0.3. 
Two distinct branches of ESs are displayed. Along each branch, the profile changes continuously. Figure HJ1 
displays an array of profiles of these ESs for D = 100 and D — 5. Note that the second (higher-(5) branch 
can be considered to be a family of fundamental, i.e. single-humped, solitons throughout the range of existence, 
whereas the double-humped structure of the SH component of the first branch becomes evident for small q. This 
property, and the location of the branches in the (6, q) plane, are fully consistent with results obtained in the 
continuum limit 0, Fig. 3]. Our numerical results also suggest that there may exist further branches for lower 
6- values, each subsequent one containing more oscillations in core of the SH component. However, computation 
becomes unreliable beyond the first two branches for the reasons given at the end of the last section, and so we 
do not display those results here. Also experience from the continuum model suggests that non-fundamental 
ESs never have a chance to be stable. 

Figure 13b) depicts continuation in D of the solutions on these two branches for q — 5. In this case, we find 
that there is a lower limit on D beyond which no ESs exist. This corresponds to the right-hand inequality in 
the first condition in (|23|) . 

We will now check the validity of the analytical approximation developed in Section 3 when k and hence k 
are large. In Fig.0 we display the two numerically computed ES branches with the same values of 71 and 72 as 
those in Fig.JSJfor D = 10 when k is rather large, and compare them with the analytical prediction (|17fl . taking 
into regard the relations (|22|l . In Fig.EJa), the parameters 8 and k are varied with q = 2k, and in Fig.[7Jb) the 
parameters 5 and q are varied for k = 50. The analytical predictions are plotted as dashed lines. We find good 
agreement between the numerical and analytical results, especially for large k. Figure |S1 displays the profiles 
of these ESs. The fundamental soliton in Fig. Eta) has a steep peak at n = 0, as assumed in the approximate 
analysis of Section 3. 

It is well known that in the continuum case multi-pulse homoclinic orbits exist (under a mild Birkhoff 
signature condition 0, |44| ) along families of curves in the parameter plane that accumulate on the curves 
of the fundamental ES solutions. Unlike the branches computed above, they do not feature a single-humped 
shaped in either component, but rather look like bound-states of two spatially separated fundamental solitons. 
We have found exactly the same solution families in the discrete model too. An example is presented in Fig. [5] 
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FIG. 5: Branches of the discrete ES solutions for k = 0.3 and 71 = 72 = 0.05. (a) Solutions in the (8, g)-plane for 
(D,N) = (100,86) (solid line) (10,31) (dashed line) and (5,23) (broken line). According to ESs also exist only 
above the dotted line q = 0.6 for this fc-value. Labeled points correspond to the relevant subpanels of Fig. [(J where 
solutions profiles are displayed, (b) Solutions in the (D, 5)-plane for q = 5 and N = 22. Note that, according to 12311 . we 
must have SD > 2.2 (the boundary is shown as a dotted line). 



B. The case of 5 < 

Figure shows ES branches in the case of self-defocusing x^ 3 ' terms, with 71 = 72 = —0.05 and S < 0. 
Figures ITTI and IT21 display the profiles of these ESs for D = 100, D = 10 and D = 6. Note that these curves 
and the solutions on them for D = 100 arc identical, to the accuracy depicted, to the corresponding ones found 
in the continuum counterpart of the model in Ref. 0, Figs. 2,3]. Notice the variety of multi-humped shapes 
of the solitons belonging to these families. Like the continuum model, only the inner-most of these branches 
represents a fundamental soliton. Continuation towards the anti-continuum limit, D — > 0, becomes numerically 
problematic for these solitons. It was found difficult to compute the solutions with repeatable accuracy while 
varying TV below D w 6. Clearly, these branches become increasingly spiky as D is decreased (see Fig. 1121 . and 
it may happen that branches of ES solutions actually terminate before they reach the minimum value of D at 
which they remain embedded, which would be D min = 2k — 1, for the values of q and S used in Fig. HOf b). 

Figure IT31 shows the ES branches in a still more exotic case of opposite signs in front of the FF and SH x^ 
SPM terms, 71 = 0.05 and 72 = —0.05. This case may, in principle, also be physically relevant - not to ordinary 
optical materials, but rather to photonic crystals (see, e.g., 0] and references therein). Figure PHI displays the 
profiles of these ESs for D = 100 and D = 10. Note similarity with the branches in Fig.^]for small 5. However, 
in this case it is found that the solution branches still exist for large k, rather than undergoing turning back with 
the increase of k. Also the internal oscillations on the non-fundamental branches become far less pronounced. 

The approximate analysis of Section 3 is also valid for 71 > and 72, 5 < when k is large. Figure IT51 shows 
the fundamental ES branch (the left one in Fig. ll3f a)) with the same values of 71 and 72 as those in Fig. 1131 for 
D = 10 when k is rather large. In Fig. I15f a). the parameters S and k are varied for q = 2k, and in Fig. I15f b) 
the parameters S and q are varied for k = 50, as in Figs. E£a) and (b). The predictions based on Eqs. I|17f) 
with ill'L'l) are plotted as dashed lines. A profile of the ES is also displayed in Fig. I15f a). We see that Eq. I|17|) 
again approximates well the numerical result for the fundamental solitons when k is large. The ES in Fig. llSf a) 
features a steep peak at n — 0, as assumed in the analytical approximation. 

Finally, we briefly discuss the case in which 71 = 72 — 0, i.e., the v/ 3 ) terms are absent. As shown in 
Fig. 1161 we could follow solutions of the three-dimensional algebraic problem for G e . However, these solutions 
do not represent ESs since the origin is not a saddle-center but a hyperbolic saddle. Although the region 
q/2 < k < q/2 — 5D in which the origin is a saddle-center becomes wider when D is larger, for the solution k 
diverges to 00 as D — * 00. Thus, there seems to be no hope that an ES exists in this case. As mentioned above, 
this situation is very similar to that in the continuum model |JT]|. 
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FIG. 6: Profiles of the discrete ESs corresponding to the labeled points in Fig. EJa): (a) 8 — 0.58637 and q = 5; (b) 
S = 2.9894 and q = 5; (c) 6 = 0.65467 and q = 0.6; (d) 6 = 3.3693 and q = 0.6; (e) 5 = 1.3496 and q = 5; (f) 6 = 3.1432 
and g = 5; (g) 5 = 0.80213 and g = 0.6; (h) S = 3.3824 and q = 0.6. In panels (a)-(d), D = 100 and N = 85, and in 
panels (e)-(h), D = 5 and iV = 22. In this and all subsequent plots, the FF (it-component) is interpolated by a solid line 
and the SH (u-component) by a broken line. 
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FIG. 7: Branches of discrete ESs for 71 = 72 = 0.05 and D = 10; N = 7 is chosen for larger 8, and N = 5 for smaller 
5. The analytical prediction given by Eqs. I|17^ and 1221 is plotted as a dashed line, (a) Solutions in the (S, fc)-plane for 
q — 2k. (b) Solutions in the (S, q)-plane for k = 50. The ESs exist in the region 100 < q < 100 + 105, whose boundary is 
shown as a dotted line. 
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FIG. 8: Examples of the ESs on solution branches in Fig.EJb): (a) 5 = 0.34718, q = 102 and N = 7; (b) 5 = 0.0029014, 
q = 100.02 and N = 5. 




FIG. 9: A typical example of a two-pulse bound state ES in the lattice system I12H with S = 1.1261, 71 = 72 = 0.05, 
D = 5, k = 0.3, 3 = 7 and N = 51. 
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FIG. 10: Branches of ESs in the discrete system @ with q — 1 and 71 = 72 = —0.05: (a) In the (S, fc)-plane for D — 100 
and 10; (b) in the (D, fc)-plane for S — —0.5 and N = 30. In panel (a), the solid and dashed curves represent the results 
for D — 100 and 10, respectively. When D — 100 (resp. D = 10), N = 85 (resp. N = 30) were used for the most part 
but ./V = 93 or N — 109 (resp. N = 35 or N = 40) were used for the two outer curves with large According to 
condition 12311 . ESs exist only above k = 0.5 and below k = 0.5 — SD, when q — 1. 



VII. CONCLUSIONS 



In this paper, we have established the existence of embedded solitons (ESs) in the discrete model of the 
second-harmonic generation in the presence of cubic nonlinearity. The model can be naturally realized as an 
array of channel waveguides in the spatial domain, therefore our results suggest possibilities for new experiments 
with discrete spatial solitons in nonlinear optics. We have also introduced a simplified model, with the x' 3 ' 
nonlinearity present solely at the central site, in which the existence of an ES was demonstrated in an asymptotic 
analytical form. In the general case, the existence region for the ESs in the discrete model was found numerically. 
Moreover, we have checked that the asymptotic analysis of the simplified model in the limit of large wavenumbers 
gives a good approximation of the codimension-one set in parameter space, on which the ESs exist. 

More generally, we have established, that unlike the discretizations of other (dissipative) continuum models 
bearing localized solutions, the codimension of ESs does not change when one passes to a discrete version. 
Whereas continuum ESs may be regarded as homoclinic orbits to saddle-center equilibria in reversible flows 
(ordinary differential equations), discrete ESs should be thought of as homoclinic orbits to saddle-center fixed 
points of reversible maps. Both have codimension one in the parameter space. Understanding this property 
led us to the derivation of a numerical method for computing homoclinic orbits to nonhyperbolic equilibria of 
reversible maps. 

Accurate investigation of stability of ESs in the lattice model is beyond the scope of the present investigation. 
Systematic results for the stability will be presented elsewhere; however, some preliminary results suggest that 
the fundamental discrete ESs inherit the semi- stability property found in the corresponding continuous model 
|HlH|. We expect that general arguments in favor of the semi-stable character of these solitons for the continuous 
case should apply in the discrete setting as well. 

Figure IT71 shows a preliminary computational result for D = 10, 6 = —1, k = 0.6895, q = 1 and 71 = 72 = 
—0.05, corresponding to the fundamental ES of Fig. Illf e). Here Eq. J2J was integrated numerically by the 
fourth-order Runge-Kutta method under the boundary condition 

u_ N _ 1 (z) = u N+1 {z) = v_ N _ 1 (z) = v N+1 (z) = (37) 

for all z > with N = 93 and the initial condition 

u n (0) = (1 + a)U n , « n (0) = (1 + c 2 )F n (38) 

where (U n ,V n ) represents an approximate ES given by the data of Fig. Illf e) for n| < N — 30 and by 
U n — V n — for |n| > N, and Ci ( 2 are small amplitudes of the initial disturbance, which were chosen to 
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FIG. 11: ESs on the branches in Fig. IToTal for S = -1: (a) k = 0.69564; (b) k = 0.71312; (c) k = 0.71383; (d) 
k = 0.71387; (e) k = 0.6895; (f) k = 0.70811; (g) k = 0.70898; (h) k = 0.70902. In panels (a)-(d), D = 100 and N = 85, 
and in panels (e)-(h), D = 10 and TV = 30. 
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FIG. 12: ESs at the end of the above two branches in Fig.EJb) for D = 6: (a) k = 0.78427; (b) k = 0.78426. 
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FIG. 13: Branches of ESs with q = 1, ji = 0.05 and 72 = -0.05: (a) In the (S, fc)-plane for D = 100 and 10; (b) in the 
(D, fc)-plane for k = 1 and N = 24. In panel (a), the solid and dashed curves represent the results for D = 100 and 10, 
respectively. When D = 100 (resp. D — 10), N — 85 (resp. N = 30) was used. In panel (b), the forth solid line from the 
bottom still exists although it is very short and almost invisible. The ESs exist only in the region 0.5 < k < 0.5 — SD 
when q = 1, whose boundary is indicated by a dotted line. 



be c\ — C2 — 0.01. The positive values of mean that the norm 

N 

E= (K| 2 + 2M 2 ) (39) 
n=-N 

of the perturbed state (in the temporal-domain optical model, it plays the role of energy) is greater than that 
of the unperturbed ES. In Fig.El we see the characteristic hallmark of semi-stability for the ES: The shapes of 
the perturbed wave are almost unchanged for a long period of t in Figs. ll7f a) and(b), and the amplitudes |tto(£)| 
and |uo(t)| exhibit small oscillations near the unperturbed values in Figs. lT7T c) and(d). However, the perturbed 
ES was destroyed in the simulations when different signs of c\ and C2 were chosen. Further investigation of 
stability and bifurcation will be the subject of future work. 

Finally, our results so far pertain only to those discrete ESs that are symmetric with respect to the involution 
R defined in Eq. (|28|l : solitons with this symmetry have an on-site central peak. It is also possible to apply 
techniques developed in this work to solutions symmetric with respect to involution R, see Eq. (|29|l , that should 
feature an inter-site central peak. In other physical context, such waves are less likely to be stable than waves 
that are centered on a lattice site j2(|. However, that understanding typically applies to regular (gap) discrete 
solitons and need not necessarily apply to ESs. We shall address this issue in future work. 
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FIG. 14: Examples of ESs on the branches in Fig. Eta) for k = 1: (a) 5 = -0.23601; (b) <5 = -0.099245; (c) 
<5 = -0.053237; (d) S = -0.033776; (e) 6 = -0.24838; (f) 5 = -0.13304; (g) 5 = -0.096003; (h) S = -0.05109. In panels 
(a)-(d), D = 100 and N = 85; in panels (e)-(h), D = 10 and N = 30. 
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FIG. 15: Branches of discrete ESs for 71 = 0.05, 72 = —0.05, D = 10 and N — 7. The analytical prediction, given by 
Eqs. I|17[l and \Y2'2t is plotted by dashed lines, (a) Solutions in the (8, fc)-plane for q — 2k. (b) Solutions in the (S, q)-plane 
for k = 50. The ESs exist in the region 100 + 105 < q < 100, whose boundary is shown by the dotted line. An example 
of the ES, for 8 = -0.31873 and q = 98, is plotted in panel (b). 
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FIG. 16: Branches of solutions of the three-dimensional algebraic problem for G e in the absence of the x nonlinearity, 
with 71 = 72 = 0, and q = 1. (a) In the (5, fc)-plane for D = 10; (b) in the (D, fc)-plane for 5 = —0.1 and iV = 15. 
In panel (a), N = 15 and 25 were used for k > 1 and k < 1, respectively. The ESs would exist only in the region 
0.5 < k < 0.5 — SD, the boundaries of which are plotted as dotted lines, when q — 1. An example of a non- embedded 
(gap) soliton, found in this case for 8 = —0.1 and k = 1.9027, is plotted in panel (a). 
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FIG. 17: Evolution of the fundamental ES (corresponding to Fig. Illf eO initiated by the nondestructive perturbation 
with ci = C2 = 0.01 in the initial conditions (1381 . The dashed lines in panels (c) and (d) are the amplitudes of the uo- 
and uo-components of the unperturbed soliton. 



References 



[1] Yang J, Malomed B A and Kaup D J 1999 Embedded solitons in second-harmonic-generating systems Phys. Rev. 
Lett. 83 1958-1961 

[2] Champneys A R, Malomed B A, Yang J and Kaup D J 2001 Embedded solitons: Solitary waves in resonance with 

the linear spectrum Physica D 152-153 340-354 
[3] Yang J, Malomed B A, Kaup D J and Champneys A R 2001 Embedded solitons: A new type of solitary wave Math. 

Comp. Sim 56 585-600 

[4] Champneys A R and Groves M D 1997 A global investigation of solitary waves solutions to a two-parameter model 

for water waves J. Fluid. Mech. 342 299-229 
[5] Fujioka J and Espinosa A 1997 Soliton-like solution of an extended NLS equation existing in resonance with linear 

dispersive waves J. Phys. Soc. Jpn. 66 2601-2607 
[6] Champneys A R, Malomed B A and M. J. Friedman 1998 Thirring solitons in the presence of dispersion Phys. Rev. 

Lett. 80 4168-4171 

[7] Atai J and Malomed B A 2001 Solitary waves in systems with separated Bragg grating and nonlinearity Phys. Rev. 
E 64 066617 

[8] Kolossovski K, Champneys A R, Buryak A and R. A. Sammut 2002 Multi-pulse embedded solitons as bound states 

of quasi-solitons Physica D 171 153-177 
[9] Espinosa-Ceron A, Fujioka J and Gomez-Rodriguez A 2003 Embedded solitons: Four-frequency radiation, front 

propagation and radiation inhibition Physica Scripta 67 314-324 
[10] Yang J K 2003 Stable embedded solitons Phys. Rev. Lett. 91 143903 

[11] Mak W C K, Malomed B A and Chu P L 2004 Symmetric and asymmetric solitons in linearly coupled Bragg gratings 
Phys. Rev. E 69 066610 

[12] Merhasin I M and Malomed B A 2004 Four-wave solitons in Bragg cross-gratings J. Optics B: Quant. Semiclass. 
Opt. 6 S323-S332 

[13] Pelinovsky D E and Yang J 2002 A normal form for nonlinear resonance of embedded solitons, Proc. Roy. Soc. Lond. 
A 458 1469-1497 

[14] Champneys A R and Kivshar Yu S 2000 Origin of multikinks in dispersive nonlinear systems Phys. Rev. E 61 
2551-2554 



20 



[15] Flach S, Zolotaryuk Y and Kladko K 1999 Moving lattice kinks and pulses: An inverse method Phys. Rev E 59 
6105-6115 

[16] Morgante A M, Johansson M, Kopidakis G and Aubry S 2002 Standing wave instabilities in a chain of nonlinear 

coupled oscillators Physica D 162 53-94 
[17] Peyrard M and Kruskal M D 1984 Kink dynamics in the highly discrete sine-Gordon system Physica D 14 88-102 
[18] Braun O M and Kivshar Yu S 1998 Nonlinear dynamics of the Frenkel-Kontorova model Phys. Rep. 306 1-108 
[19] Savin A V. Zolotaryuk Y and Eilbeck J C 2000 Moving kinks and nanopterons in the nonlinear Klein-Gordon lattice 

Physica D 138 267-281 

[20] Aigner A A, Champneys A R and Rothos V R 2003 A new barrier to the existence of moving kinks in Frenkel- 
Kontorova lattices Physica D 186 148-170 

[21] Gonzalez-Perez-Sandi S, Fujioka J and Malomed B A 2004 Embedded solitons in dynamical lattices Physica D 197 
86-100 

[22] Fiedler B and Scheurle J 1996 Discretization of Homoclinic Orbits, Rapid Forcing and "Invisible" Chaos Memoirs 

of Amer. Math. Soc. 119 (Providence RI: American Mathematical Society) 
[23] Bang O, Clausen C B, Christiansen P L and Torner L 1999 Engineering competing nonlinearities Opt. Lett. 24 

1413-1415 

[24] Corney J F and Bang O 2001 Solitons in quadratic nonlinear photonic crystals Phys. Rev. E 64 047601 
[25] Johansen S K, Carrasco S, Torner L and Bang O 2002 Engineering of spatial solitons in two-period QPM structures 
Opt. Commun. 203 393-402 

[26] Etrich C, Lederer F, Malomed B A, Peschel T, and Peschel U 2000 Optical solitons in media with a quadratic 

nonlinearity Progr. Opt. 41 483-568 
[27] Iwanow R, Schiek R, Stegeman G I, Pertsch T, Lederer F, Min Y and Sohler W 2004 Observation of discrete 

quadratic solitons Phys. Rev. Lett. 93 113902 
[28] Kcvrekidis P G Malomed B A, and Musslimani Z 2003 Discrete gap solitons in a diffraction-managed waveguide 

array Eur. Phys. J. D 67 013605 
[29] Eisenberg H S, Silberberg Y, Morandotti R, Boyd A and Aitchison J S 2000 Diffraction management Phys. Rev. 

Lett. 85 1863-1866 

[30] Pertsch T, Zentgraf T, Peschel U and Lederer F 2002 Anomalous refraction and diffraction in discrete optical systems 

Phys. Rev. Lett. 88 093901 
[31] Devaney R L 1976 Reversible diffeomorphisms and flows Trans. Amer. Math. Soc. 218 89-113 

[32] Lamb J S W and Roberts JAG 1998 Time-reversal symmetry in dynamical systems: A survey Physica D 112 1-39 
[33] Guckenheimer J and Holmes P 1983 Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields 
(New York: Springer) 

[34] Doedel E, Champneys A R, Fairgrieve T F, Kuznetsov Y A, Sandstede B and Wang X 1997 AUT097: Continuation 
and Bifurcation Software for Ordinary Differential Equations (with HomCont) Available by anonymous ftp from 
ftp.cs.concordia.ca, directory pub/doedel/auto 
[35] Nusse E H and Yorke J A 1997 Dynamics: Numerical Explorations 2nd ed. (New York: Springer) 
[36] Kawakami H 1981 Algorithme numerique definissant la bifurcation d'un point homocline C. R. Acad. Sc. Paris, 
Seme I 293 401-403 

[37] Beyn W-J and Kleinkauf J N 1997 The numerical computation of homoclinic orbits for maps SIAM J. Num. Anal. 
34 1207-1236 

[38] Beyn W-J and Kleinkauf J N 1997 Numerical approximation of homoclinic chaos Numer. Algorithms 14 25-53 
[39] Yagasaki K 1998 Numerical detection and continuation of homoclinic points and their bifurcations for maps and 

periodically forced systems Int. J. Bifurcation Chaos 7 1617-1627 
[40] Bergamin J M, Bountis T and Vrahatis M N 2002 Homoclinic orbits of invertible maps Nonlinearity 15 1603-1619 
[41] Yagasaki K 1998 HomMap: An Auto driver for homoclinic bifurcation analysis of maps and periodically forced 

systems (Gifu Japan: Gifu University) 
[42] Lombardi E 2000 Oscillatory Integrals and Phenomena beyond All Algebraic Orders: With Applications to Homoclinic 

Orbits in Reversible Systems (Berlin: Springer) 
[43] Champneys A R 2001 Codimension-one persistence beyond all orders of homoclinic orbits to singular saddle centers 

in reversible systems Nonlinearity 14 87-112 
[44] Mielke A, Holmes P and O'Reilly O 1992 Cascades of homoclinic orbits to, and chaos near, a Hamiltonian saddle- 
center J. Dyn. Diff. Eqns. 4 95-126 
[45] Andreani L C, Agio M, Bajoni D, Belotti M, Galli M, Guizzetti G, Malvezzi A M, Marabelli F, Patrini M and 

Vecchi G 2003 Morphology and optical properties of bare and polydiacetylenes-infiltrated opals Synthetic Metals 

139 695-700 



21 



